clc
clear

global NUMBER_OF_EVENTS LENGTH BIG TNOW MAINT_COST REPAIR_COST SALES_REVENUE EXP_RATE 
global SHAPE1 SCALE1 SHAPE2 SCALE2 A3 B3 B U L MU STATES_MAX THRESHOLD


init_global1




for rep=1:8
% initialize event clock and the structure 
rng(22+rep*10,'twister');
for i=1:NUMBER_OF_EVENTS
    event_clock(i)=0;
end


% revenue contains the actual revenue since last update and time contains
% the actual time when last update was performed
% TR and TT contain the total reward and total time needed by the algorithm
% TR and TT are not to be confused with actual time since they are updated only when a non-exploratory action occurs 

alg=struct('iter',0,'old_action',1,'prodcount',1,'old_prodcount',0,'revenue',0,'time',0,...
'time_for_failure',gamrnd(SHAPE1,SCALE1),...
'rho',0,'flag',1,'TR',0,'TT',0,'Pi',zeros(STATES_MAX,2),'J',zeros(STATES_MAX));


% initialize values to event clock

event_clock(1)=LENGTH;
event_clock(2)=12; % time of next production 
event_clock(3)=BIG; % time of next repair: unknown 
event_clock(4)=BIG; % time of next maintenance: unknown

% simulation begins


done = 0; % Pnemonic for simulation, 1 stands for end 0 stands for continue 

    while 0 == done
                
        next_event=time_setter(event_clock);
        
        % Schedule the next event  
			
        switch next_event
                
                case 1
                % end of simulation
                done = 1;
                
	            case 2
                % produce
                [event_clock,alg]=produce(event_clock,alg);
                
                case 3
                % repair
                [event_clock,alg]=repair(event_clock,alg);
                
                case 4
                % maintain
                [event_clock,alg]=maintain(event_clock,alg);
                
				
                
                
        end % End of switch statements */	

   end % of while loop
                 
 rho=alg.TR/TNOW

 avg_reward(rep)=rho;
 end
 tvalue=2.36;
 avg=mean(avg_reward)
 std_dev=std(avg_reward);
 hw=tvalue*std_dev/sqrt(8)
 CI_U=avg+hw
 CI_L=avg-hw
   
 
		
